#include <cv.h>
#include <highgui.h>
#include <math.h>
#include "fftw3.h"

#define N 300

using namespace cv;

int round_gint (double value)
{
  double floored = floor (value);
  if (value - floored > 0.5)
  {
    return (int)(floored + 1);
  }
  return (int)floored;
}

int boost (double value)
{
  double bounded = fabs (value / 160.0);
  int boosted = round_gint (128.0 * sqrt (bounded));
  boosted = (value > 0) ? boosted : -boosted;
  return boosted;
}

int reorder_row (int row, int height)
{
  int h2 = (height + 1)/2;
  int row2 = (row % 2 == 0) ? h2 - row/2 - 1 : h2 + row/2;
  row2 = (row2 + height) % height;
  return row2;
}


int reorder_col (int col, int width)
{
  return (col + (width + 1)/2 - 1) % width;
}

double normalize (int x, int y, int width, int height)
{
  double cx = (double)abs(x - (width + 1)/2 + 1);
  double cy = (double)abs(y - (height + 1)/2 + 1);
  double energy = (sqrt(cx)+sqrt(cy));
  return energy*energy;
}

int fft(Mat *ref, Mat *res)
{
	int 	i, j, k;
	int col2,row2,bounded;
	double	tmp,v,norm;
	CV_Assert(ref->depth() != sizeof(uchar)); 
	/* get image properties */
	int width    = ref->cols;
	int height   = ref->rows;
	int step     = ref->step;
	int fft_size = width * height;
        printf("%u %u %u\n",width,height,step);


	/* setup pointers to images */
	uchar 	*ref_data = ( uchar* ) ref->data;
	uchar   *res_data = ( uchar* ) res->data;



	
	/* allocate FFTW input and output arrays */
	fftw_complex *img1 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * fft_size );
	fftw_complex *img2 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * fft_size );


	/* setup FFTW plans */
	//fftw_plan fft_img1 = fftw_plan_dft_1d( fft_size, img1, img1, FFTW_FORWARD,  FFTW_ESTIMATE );
	fftw_plan fft_img1 = fftw_plan_dft_2d( width , height, img1, img1, FFTW_FORWARD,  FFTW_ESTIMATE );
	fftw_plan fft_img2 = fftw_plan_dft_2d( width , height, img1, img1, FFTW_BACKWARD,  FFTW_ESTIMATE );

		
	//load images' data to FFTW input 
	for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++, k++ ) {
			img1[k][0] = (double) ref_data[i * step + j];
			img1[k][1] = 0.0;			
		}
	}
	
	/* obtain the FFT of img1 */
	fftw_execute( fft_img1 );
	//fftw_execute( fft_img2 );

	//normalize and copy to result image
	for( i = 0, k = 0 ; i < height ; i++ ) {
		//printf("%d:",i );
                row2 = reorder_row(i, height);
		for( j = 0 ; j < width ; j++, k++ ) {
		        col2 = reorder_col(j, width);
		        norm = normalize(col2,row2, width, height);
  		        v = img1[k][0]/(double)fft_size;
		        bounded = boost (v*norm);
		        res_data[row2 * step + col2]=  bounded;
			//printf("%u,",res_data[i * step + j]);
			//printf("%u %u %f %f %u %u,",row2,col2,norm,v,bounded,res_data[row2 * step + col2]);
		}
		//printf("\n");
	}


        
	/*
	for( i = 0 ; i < fft_size ; i++ ) {
        	res_data[i] = (uchar) (img2[i][0]/ ( double )fft_size);
		if (img2[i][0]!=0 ) {printf("%f %u\n",img2[i][0],res_data[i]);}
		//ref_data[i] = img1[i][0] ;
	}*/

	/* deallocate FFTW arrays and plans */
	fftw_destroy_plan( fft_img1 );
	fftw_free( img1 );

}

int normalizeFFT( Mat *ref, Mat *res)
{
	int 	i, j, k;
	double	tmp;
	CV_Assert(ref->depth() != sizeof(uchar)); 
	/* get image properties */
	int width    = ref->cols;
	int height   = ref->rows;
	int step     = ref->step;

        printf("%u %u %u\n",width,height,step);
	uchar 	*ref_data = ( uchar* ) ref->data;
	uchar   *res_data = ( uchar* ) res->data;

	for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++, k++ ) {
		        tmp +=  abs(ref_data[i * step + j]);
			//printf("%d %u\n",i * step + j,res_data[i * step + j]);
		}
	}
	printf("suma:%f\n",tmp);

	for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++, k++ ) {
		        res_data[i * step + j]=  (double) abs(ref_data[i * step + j])/tmp;
			//printf("%d %u\n",i * step + j,res_data[i * step + j]);
		}
	}

}

int entropy(Mat img)
{
    //cvtColor(img, hsv, CV_BGR2HSV);
    //cvtColor(img, res, CV_BGR2HSV);
    //cvtColor( img, res, CV_RGB2GRAY );
    //cvtColor(img, hsv, CV_BGR2Lab);
    //cvtColor(img, hsv, CV_RGB2Lab);
    //cvtColor(img, hsv, CV_RGB2HSV);


    Mat res;
    res = Mat::eye(img.rows,img.cols, CV_8UC1)*255; 	
    fft(&img,&res);
    //imshow( "img", img );
    //imshow( "fft", res );
    //normalizeFFT(&res,&img);
    //imshow( "FFTN", img );
    //imwrite("kk.png",res);		
    //waitKey(0); 


    // Quantize the hue to 16 levels
    int hbins = 255;
    int histSize[] = {hbins};
    // varies from 0 to 179, see cvtColor
    float hranges[] = { 0, 255 };
    const float* ranges[] = { hranges };

    MatND hist;
    // we compute the histogram from the 0-th and 1-st channels

    int channels[] = {0};

    calcHist( &res, 1, channels, Mat(), // do not use mask
             hist, 1, histSize, ranges,
             true, // the histogram is uniform
             false );

     // Draw the histograms for R, G and B
 int hist_w = 400; int hist_h = 400;
 int bin_w = cvRound( (double) hist_w/hbins );

 Mat histImage( hist_w, hist_h, CV_8UC3, Scalar( 0,0,0) );

 /// Normalize the result to [ 0, histImage.rows ]
 normalize(hist,hist, 0, histImage.rows, NORM_MINMAX, -1, Mat() );

 /// Draw for each channel
 for( int i = 1; i < hbins; i++ )
   {
     line( histImage, Point( bin_w*(i-1), hist_h - cvRound(hist.at<float>(i-1)) ) ,
                      Point( bin_w*(i), hist_h - cvRound(hist.at<float>(i)) ),
                      Scalar( 0,255,0), 2, 8, 0  );
    }

 /// Display
 namedWindow("calcHist Demo", CV_WINDOW_AUTOSIZE );
 imshow("calcHist Demo", histImage );

// waitKey(0);


    //ENTROPY
    double H=0;	
    for( int h = 0; h < hbins; h++ ) 
           {
            float binVal = hist.at<float>(h);
	    if (binVal != 0) {
	   	    H+=binVal*log(float(binVal));	
	    }
        }
    printf("ENTROPY:\%f\n",H);

    return H;
}


int main( int argc, char** argv )
{
    Mat src;	
    double H[N];
    double min=1000000000000,max=0;
    for(int i=0;i<N;i++)
    {
        char filename[150];
        sprintf(filename, "../img/luna-%d.png", (i+1));
        //sprintf(filename, "../img/polio%d.JPG", (i));
	//sprintf(filename, "../img/sin2.png");
        //if(!(src=imread(filename, 1)).data ) return -1; //COLOR
        if(!(src=imread(filename, 0)).data ) return -1; //GREYSCALE
        H[i]=entropy(src);
        //printf("%d :  %f\n",i+1,H[i]);
        if (min >=H[i]) min=H[i]; 
        if (max <=H[i]) {
		max=H[i];
                printf("%d :  %f\n",i+1,max);
	}
    }
    //NORMALIZE to 100
    double h,minmax=max-min;
    for(int i=0;i<N;i++)
    {
       h=100-(H[i]-min)*100/minmax;
       h=H[i];
       printf("%d :  %f\n",i+1,h);
    }
	

/*    if( argc != 2 || !(src=imread(argv[1], 1)).data )
        return -1; */

}



